################################################################
### This file contains code to do the synthetic control
### analysis in the paper "The Hirschman Effect of International
### Financial Ties", Xun Pang and Chong Chen, 2020
################################################################
#setwd("Replication")
rm(list=ls())
library(lme4)
library(arm)
library(plm)
library(dplyr)
library(texreg)
library(extrafont)
library(ggplot2)
library(showtext)
library(ggthemes)
library(scales)
font_add('KaiTi', 'KaiTi.ttf')
showtext_auto(enable = TRUE)
######### fixed-effect model with SD clustered by ccode

source("Code/setup.R")

# Note: all Dvs have been led for one period (which is equivalent to lagging all IVs )
load("data.RData")


## Let's call the models mixed-effect models (economists like it more)
summary(data$year)
data <- data %>% 
        dplyr::filter(ccode !=710 )

data <- data %>% 
  dplyr::mutate(defense_usa_lag1 = ifelse(is.na(defense_usa_lag1), 0, defense_usa_lag1)) %>% 
  dplyr::filter(year < 2019 & year > 2008) 


#distance_idealpoint_chnlog
#tradedepp: trade dependence on the world

## MLM

# m0: simple model
m0 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1','swap_dummy_lag1','distance_idealpoint_usalog_lag1','minidist_usalog_lag1',
                       'defense_usa_lag1'),
               modname = c('base'),
               data = data)

summary(m0$base)


m1 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('swap_dummy_lag1',
                       'polity2_lag1', 'gdppc_log_lag1', 'pop_log_lag1','milexp_lag1',
                       'sco_full_lag1','asean_memb_lag1'),
               modname = c('country'),
               data = data)

summary(m1$country)


m2 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1', 'swap_dummy_lag1',
                       'dep_export_lag1','dep_import_lag1', 
                       'arms_exportfromCHNlog_lag1', 'arms_importtoCHNlog_lag1',
                       'bits_chn_lag1','pta_china_lag1',
                       'minidist_log_lag1','partner_level_lag1'),
               modname = c('bilateral'),
               data = data)

summary(m2$bilateral)

m3 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1', 'swap_dummy_lag1',
                       'polity2_lag1', 'gdppc_log_lag1', 'pop_log_lag1','milexp_lag1',
                       'dep_export_lag1','dep_import_lag1', 
                       'arms_exportfromCHNlog_lag1', 'arms_importtoCHNlog_lag1',
                       'bits_chn_lag1','pta_china_lag1','sco_full_lag1','asean_memb_lag1',
                       'minidist_log_lag1','partner_level_lag1'),
               modname = c('full'),
               data = data)

summary(m3$full)

m4 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1', 'swap_dummy_lag1',
                       'polity2_lag1', 'gdppc_log_lag1', 'pop_log_lag1','milexp_lag1',
                       'dep_export_lag1','dep_import_lag1', 
                       'arms_exportfromCHNlog_lag1', 'arms_importtoCHNlog_lag1',
                       'bits_chn_lag1','pta_china_lag1','sco_full_lag1','asean_memb_lag1',
                       'minidist_log_lag1','partner_level_lag1',
                       'distance_idealpoint_usalog_lag1','minidist_usalog_lag1',
                       'defense_usa_lag1'),
               modname = c('us'),
               data = data)

summary(m4$us)

###############

m0_1 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1','swap_rmbsizelog_lag1','distance_idealpoint_usalog_lag1','minidist_usalog_lag1',
                       'defense_usa_lag1'),
               modname = c('base'),
               data = data)

summary(m0_1$base)

m1_1 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('swap_rmbsizelog_lag1',
                       'polity2_lag1', 'gdppc_log_lag1', 'pop_log_lag1','milexp_lag1',
                       'sco_full_lag1','asean_memb_lag1'),
               modname = c('country'),
               data = data)

summary(m1$country)



m2_1 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1', 'swap_rmbsizelog_lag1',
                       'dep_export_lag1','dep_import_lag1', 
                       'arms_exportfromCHNlog_lag1', 'arms_importtoCHNlog_lag1',
                       'bits_chn_lag1','pta_china_lag1',
                       'minidist_log_lag1','partner_level_lag1'),
               modname = c('bilateral'),
               data = data)

summary(m2_1$bilateral)

m3_1 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1', 'swap_rmbsizelog_lag1',
                       'polity2_lag1', 'gdppc_log_lag1', 'pop_log_lag1','milexp_lag1',
                       'dep_export_lag1','dep_import_lag1', 
                       'arms_exportfromCHNlog_lag1', 'arms_importtoCHNlog_lag1',
                       'bits_chn_lag1','pta_china_lag1','sco_full_lag1','asean_memb_lag1',
                       'minidist_log_lag1','partner_level_lag1'),
               modname = c('full'),
               data = data)

summary(m3_1$full)

m4_1 <- MLMs_fit(dvs = c('distance_idealpoint_chnlog'),
               ivs = c('distance_idealpoint_chnlog_lag1', 'swap_rmbsizelog_lag1',
                       'polity2_lag1', 'gdppc_log_lag1', 'pop_log_lag1','milexp_lag1',
                       'dep_export_lag1','dep_import_lag1', 
                       'arms_exportfromCHNlog_lag1', 'arms_importtoCHNlog_lag1',
                       'bits_chn_lag1','pta_china_lag1','sco_full_lag1','asean_memb_lag1',
                       'minidist_log_lag1','partner_level_lag1',
                       'distance_idealpoint_usalog_lag1','minidist_usalog_lag1',
                       'defense_usa_lag1'),
               modname = c('us'),
               data = data)

summary(m4_1$us)


##### Appendix Table 2-3


texreg(list(m0$base, m1$country, m2$bilateral,  m3$full,  m4$us),
       file = "./tables/tabdummy.tex",
       stars = c(0.01, 0.05, 0.1),
       digits = 3,
       scalebox = .7,
       dcolumn = FALSE, use.packages = TRUE, longtable = FALSE,
       caption = "Mixed-effects model estimations for the effects of swaps on foreign policy convergence (swap dummy)", caption.above = TRUE,
       label = "tab:swap")

htmlreg(list(m0$base, m1$country, m2$bilateral,  m3$full,  m4$us), 
        file = "./tables/tabdummy.doc", 
        stars = c(0.01, 0.05, 0.1),
        digits = 3,
        dcolumn = FALSE, use.packages = TRUE, longtable = FALSE,
        inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE)


texreg(list(m0_1$base, m1_1$country, m2_1$bilateral, m3_1$full, m4_1$us),
       file = "./tables/tabsize.tex",
       stars = c(0.01, 0.05, 0.1),
       digits = 3,
       scalebox = .7,
       dcolumn = FALSE, use.packages = TRUE, longtable = FALSE,
       caption = "Mixed-effects model estimations for the effects of swaps on foreign policy convergence (swap size)", caption.above = TRUE,
       label = "tab:swapsize")

### need to close word app
htmlreg(list(m0_1$base, m1_1$country, m2_1$bilateral, m3_1$full, m4_1$us), 
        file = "./tables/tabsize.doc", 
        stars = c(0.01, 0.05, 0.1),
        digits = 3,
        dcolumn = FALSE, use.packages = TRUE, longtable = FALSE,
        inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE)


save(m0, m1, m2, m3, m4, m0_1, m1_1, m2_1, m3_1, m4_1, file = "results.RData")



################ ################################################Figure 3################################

###################  load the results

load("results.RData")


fig_marginal_m1 <-  MLM_maginal_plot(Fits = list(m0$base, m1$country, m2$bilateral,  m3$full,  m4$us), n_sim = 1000,
                                     modname = c('m0','m1','m2','m3','m4'),
                                var = c('swap_dummy_lag1',
                                        'swap_dummy_lag1',
                                        'swap_dummy_lag1',
                                        'swap_dummy_lag1',
                                        'swap_dummy_lag1'),
                                val1 = 0, val2 = 1,
                                data = data)

fig_marginal_m1$plot +  
          labs(title = "",
               x = "平均边际效应") + 
        scale_y_discrete(labels = c("模型 5", "模型 4", "模型 3", "模型 2", "模型 1"))+
        theme(text = element_text(size=12, family ='STSong'))
  
ggsave("graph/marginal_swap.png", width = 7.5, height = 5, dpi = 500)


fig_marginal_m2 <-  MLM_maginal_plot(Fits = list(m0_1$base, m1_1$country, m2_1$bilateral, m3_1$full, m4_1$us), n_sim = 1000,
                                     modname = c('m0','m1','m2','m3','m4'),
                                     var = c('swap_rmbsizelog_lag1',
                                             'swap_rmbsizelog_lag1',
                                             'swap_rmbsizelog_lag1',
                                             'swap_rmbsizelog_lag1',
                                             'swap_rmbsizelog_lag1'),
                                     val1 = 0, val2 = 1,
                                     data = data)

fig_marginal_m2$plot +  
  labs(title = "",
       x = "平均边际效应") + 
  scale_y_discrete(labels = c("模型 5", "模型 4", "模型 3", "模型 2", "模型 1"))+
  theme(text = element_text(size=12, family ='STSong'))
ggsave("graph/marginal_swapsize.png", width = 7.5, height = 5,dpi = 500)

################################################################################################################
################ Figure 2#################
################################################################################################################

varnames =lapply(list(m0$base, m1$country, m2$bilateral,  m3$full,  m4$us), function(x) FUN = names(fixef(x)))
unique(unlist(varnames))

p <- MLM_coef(modResults = list(m0$base, m1$country, m2$bilateral,  m3$full,  m4$us), data = data, 
            vars = c("(Intercept)", 
                     "minidist_usalog_lag1",
                     "defense_usa_lag1",
                     "polity2_lag1",
                     "gdppc_log_lag1",
                     "pop_log_lag1",
                     "milexp_lag1" ,
                     "sco_full_lag1",
                     "asean_memb_lag1",
                     "dep_export_lag1",
                     "dep_import_lag1",
                     "arms_exportfromCHNlog_lag1",
                     "arms_importtoCHNlog_lag1",
                     "bits_chn_lag1",
                     "pta_china_lag1",
                     "minidist_log_lag1", 
                     "partner_level_lag1", 
                     "distance_idealpoint_usalog_lag1",
                     "distance_idealpoint_chnlog_lag1",
                     "swap_dummy_lag1")) + 
  scale_colour_manual(values = c("gold", "red", "pink", "black", "blue"))  + 
    ggtitle("混合效应模型回归系数图：货币互换协定")



p  + scale_x_discrete(labels = parse(text = c("截距", 
                                 "与美国的地理距离[list(log,t-1)]",
                                 "与美国防御性盟约[list(t-1)]",
                                 "政体分值[list(t-1)]",
                                 "人均国民生产总值[list(log,t-1)]",
                                 "人口总数[list(log,t-1)]",
                                 "军费开支比重[list(t-1)]" ,
                                 "上合组织成员国[list(t-1)]",
                                 "东盟成员国[list(t-1)]",
                                 "对中国出口贸易依赖度[list(t-1)]",
                                 "对中国进口贸易依赖度[list(t-1)]",
                                 "对中军售出口总额[list(log,t-1)]",
                                 "对中军事进口总额[list(log,t-1)]",
                                 "双边投资协议[list(t-1)]",
                                 "贸易优惠协定[list(t-1)]",
                                 "与中国的地理距离[list(log,t-1)]", 
                                 "中国外交伙伴关系类型[list(t-1)]", 
                                 "与美国的外交立场距离[list(log,t-1)]",
                                 "与中国的外交立场距离[list(log,t-1)]",
                                 "货币互换协定[list(t-1)]")))

ggsave("graph/coef1.png", width = 9, height = 10)

varnames =lapply(list(m0_1$base, m1_1$country, m2_1$bilateral, m3_1$full, m4_1$us), function(x) FUN = names(fixef(x)))
unique(unlist(varnames))

p2 <- MLM_coef(modResults = list(m0_1$base, m1_1$country, m2_1$bilateral, m3_1$full, m4_1$us), data = data, 
              vars = c("(Intercept)", 
                       "minidist_usalog_lag1",
                       "defense_usa_lag1",
                       "polity2_lag1",
                       "gdppc_log_lag1",
                       "pop_log_lag1",
                       "milexp_lag1" ,
                       "sco_full_lag1",
                       "asean_memb_lag1",
                       "dep_export_lag1",
                       "dep_import_lag1",
                       "arms_exportfromCHNlog_lag1",
                       "arms_importtoCHNlog_lag1",
                       "bits_chn_lag1",
                       "pta_china_lag1",
                       "minidist_log_lag1", 
                       "partner_level_lag1", 
                       "distance_idealpoint_usalog_lag1",
                       "distance_idealpoint_chnlog_lag1",
                       "swap_rmbsizelog_lag1")) + 
  scale_colour_manual(values = c("gold", "red", "pink", "black", "blue"))  + 
  ggtitle("混合效应模型回归系数图：货币互换协定规模")

p2 <- p2 + scale_x_discrete(labels = parse(text = c("截距", 
                                                    "与美国的地理距离[list(log,t-1)]",
                                                    "与美国防御性盟约[list(t-1)]",
                                                    "政体分值[list(t-1)]",
                                                    "人均国民生产总值[list(log,t-1)]",
                                                    "人口总数[list(log,t-1)]",
                                                    "军费开支比重[list(t-1)]" ,
                                                    "上合组织成员国[list(t-1)]",
                                                    "东盟成员国[list(t-1)]",
                                                    "对中国出口贸易依赖度[list(t-1)]",
                                                    "对中国进口贸易依赖度[list(t-1)]",
                                                    "对中军售出口总额[list(log,t-1)]",
                                                    "对中军事进口总额[list(log,t-1)]",
                                                    "双边投资协议[list(t-1)]",
                                                    "贸易优惠协定[list(t-1)]",
                                                    "与中国的地理距离[list(log,t-1)]", 
                                                    "中国外交伙伴关系类型[list(t-1)]", 
                                                    "与美国的外交立场距离[list(log,t-1)]",
                                                    "与中国的外交立场距离[list(log,t-1)]",
                                                    "货币互换协定规模[list(t-1)]")))
ggsave("graph/coef2.png", plot = p2, width = 9, height = 10)


##################### Appendix Figure 1 #################################################
###############################################################
library(ggcorrplot)

cor_df <- data[,c("minidist_usalog_lag1",
                  "defense_usa_lag1",
                  "polity2_lag1",
                  "gdppc_log_lag1",
                  "pop_log_lag1",
                  "milexp_lag1" ,
                  "sco_full_lag1",
                  "asean_memb_lag1",
                  "dep_export_lag1",
                  "dep_import_lag1",
                  "arms_exportfromCHNlog_lag1",
                  "arms_importtoCHNlog_lag1",
                  "bits_chn_lag1",
                  "pta_china_lag1",
                  "minidist_log_lag1", 
                  "partner_level_lag1", 
                  "distance_idealpoint_usalog_lag1",
                  "distance_idealpoint_chnlog_lag1",
                  "swap_dummy_lag1")]

cor_df <- na.omit(cor_df)

corr <- round(cor(cor_df), 1)

ggcorrplot(corr, hc.order = F, type = "lower",  colors = c("blue", "white", "red"),
           lab = TRUE, sig.level = 0.05, ggtheme = ggplot2::theme_minimal) + 
  theme(legend.title=element_blank(),
        axis.text = element_text(size=14),
        text = element_text(size=14,family ='STSong'),
        plot.title = element_text(hjust = .5, size = 14, face = "bold"))+
       scale_y_discrete(labels = parse(text = c( "与美国的地理距离[list(log,t-1)]",
                                                 "与美国防御性盟约[list(t-1)]",
                                                 "政体分值[list(t-1)]",
                                                 "人均国民生产总值[list(log,t-1)]",
                                                 "人口总数[list(log,t-1)]",
                                                 "军费开支比重[list(t-1)]" ,
                                                 "上合组织成员国[list(t-1)]",
                                                 "东盟成员国[list(t-1)]",
                                                 "对中国出口贸易依赖度[list(t-1)]",
                                                 "对中国进口贸易依赖度[list(t-1)]",
                                                 "对中军售出口总额[list(log,t-1)]",
                                                 "对中军事进口总额[list(log,t-1)]",
                                                 "双边投资协议[list(t-1)]",
                                                 "贸易优惠协定[list(t-1)]",
                                                 "与中国的地理距离[list(log,t-1)]", 
                                                 "中国外交伙伴关系类型[list(t-1)]", 
                                                 "与美国的外交立场距离[list(log,t-1)]",
                                                 "与中国的外交立场距离[list(log,t-1)]")))+ 
  scale_x_discrete(labels = parse(text = c( "与美国防御性盟约[list(t-1)]",
                                            "政体分值[list(t-1)]",
                                            "人均国民生产总值[list(log,t-1)]",
                                            "人口总数[list(log,t-1)]",
                                            "军费开支比重[list(t-1)]" ,
                                            "上合组织成员国[list(t-1)]",
                                            "东盟成员国[list(t-1)]",
                                            "对中国出口贸易依赖度[list(t-1)]",
                                            "对中国进口贸易依赖度[list(t-1)]",
                                            "对中军售出口总额[list(log,t-1)]",
                                            "对中军事进口总额[list(log,t-1)]",
                                            "双边投资协议[list(t-1)]",
                                            "贸易优惠协定[list(t-1)]",
                                            "与中国的地理距离[list(log,t-1)]", 
                                            "中国外交伙伴关系类型[list(t-1)]", 
                                            "与美国的外交立场距离[list(log,t-1)]",
                                            "与中国的外交立场距离[list(log,t-1)]",
                                            "货币互换协定[list(t-1)]")))

ggsave("graph/corr.png", width = 10, height = 10)



##### Appendix Table 5



library(sampleSelection)

s1 <- selection(swap_dummy ~ dep_export_lag1 + dep_import_lag1 + arms_exportfromCHNlog_lag1 + arms_importtoCHNlog_lag1 +  bits_chn_lag1 + 
               pta_china_lag1 +  sco_full_lag1 + asean_memb_lag1,
          distance_idealpoint_chnlog ~  distance_idealpoint_chnlog + distance_idealpoint_chnlog_lag1 + polity2_lag1 + gdppc_log_lag1 +
            pop_log_lag1 + milexp_lag1 +dep_export_lag1 + dep_import_lag1 + arms_exportfromCHNlog_lag1 + arms_importtoCHNlog_lag1 +
            bits_chn_lag1 + pta_china_lag1 + sco_full_lag1 + asean_memb_lag1 +
            minidist_log_lag1 + partner_level_lag1 + distance_idealpoint_usalog_lag1 + minidist_usalog_lag1 + defense_usa_lag1,
          data = data)

summary(s1)

###coefficient plot
s2 <- heckit(swap_dummy ~ dep_export_lag1 + dep_import_lag1 + arms_exportfromCHNlog_lag1 + arms_importtoCHNlog_lag1 +  bits_chn_lag1 + 
                pta_china_lag1 +  sco_full_lag1 + asean_memb_lag1,
              distance_idealpoint_chnlog ~  distance_idealpoint_chnlog + distance_idealpoint_chnlog_lag1 + polity2_lag1 + gdppc_log_lag1 +
                pop_log_lag1 + milexp_lag1 +dep_export_lag1 + dep_import_lag1 + arms_exportfromCHNlog_lag1 + arms_importtoCHNlog_lag1 +
                bits_chn_lag1 + pta_china_lag1 + sco_full_lag1 + asean_memb_lag1 +
                minidist_log_lag1 + partner_level_lag1 + distance_idealpoint_usalog_lag1 + minidist_usalog_lag1 + defense_usa_lag1,
              data = data)
summary(s2)

htmlreg(list(s1, s2), 
        file = "tables/selection2.doc", 
        stars = c(0.01, 0.05, 0.1),single.row = TRUE,
        digits = 3,
        dcolumn = FALSE, use.packages = TRUE, longtable = FALSE,
        inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE)







